SOLVE_BVP

Overview

The SOLVE_BVP function solves boundary value problems (BVPs) for systems of second-order ordinary differential equations (ODEs). Unlike initial value problems where all conditions are specified at a single point, BVPs specify conditions at two different points—typically the endpoints of an interval. This makes them essential for modeling physical systems with constraints at both ends, such as beam deflection, heat distribution in a rod, or steady-state fluid flow.

This implementation wraps SciPy’s solve_bvp function, which uses a 4th-order collocation algorithm with adaptive mesh refinement. The algorithm is based on the work of Kierzenka and Shampine, who developed the MATLAB bvp4c solver. For detailed mathematical foundations, see Numerical Solution of Boundary Value Problems for Ordinary Differential Equations by Ascher, Mattheij, and Russell.

The function solves systems of the form:

\mathbf{y}'' = A \cdot \mathbf{y}

where \mathbf{y} is a vector of dependent variables and A is a coefficient matrix. Boundary conditions are specified as \mathbf{y}(a) = \mathbf{y}_a at the left endpoint and \mathbf{y}(b) = \mathbf{y}_b at the right endpoint. Internally, the second-order system is converted to an equivalent first-order system by introducing the substitution \mathbf{u} = [\mathbf{y}; \mathbf{y}'], resulting in:

\mathbf{u}' = \begin{bmatrix} \mathbf{y}' \\ A \cdot \mathbf{y} \end{bmatrix}

The collocation method discretizes the domain and solves the resulting nonlinear system using a damped Newton iteration with an affine-invariant criterion. The solver adaptively refines the mesh until the residual tolerance (default 10^{-3}) is satisfied or the maximum number of nodes (1000) is reached. The solution is returned as a C1-continuous cubic spline evaluated at the final mesh points.

For more information on the underlying algorithm and SciPy’s implementation, see the SciPy integrate documentation and the solve_bvp source code on GitHub.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=SOLVE_BVP(y_double_prime_coeffs, bc_ya, bc_yb, x_mesh, y_init)

y_double_prime_coeffs (list[list], required): Square coefficient matrix A in the equation y’’ = A @ y (must be n×n)
bc_ya (list[list], required): Boundary condition values at the left endpoint x[0] (column vector with n rows)
bc_yb (list[list], required): Boundary condition values at the right endpoint x[-1] (column vector with n rows)
x_mesh (list[list], required): Initial mesh points for the solver (column vector, strictly increasing, minimum 2 points)
y_init (list[list], required): Initial guess for the solution (2n rows × m columns, where m is the number of mesh points)

Returns (list[list]): 2D list [x, y1, y2, …], or error message string.

Examples

Example 1: Harmonic oscillator with cosine solution

Inputs:

y_double_prime_coeffs	bc_ya	bc_yb	x_mesh	y_init
-1	1	0	0	1	0.7071	0
			0.7854	0	-0.7071	-1
			1.5708

Excel formula:

=SOLVE_BVP({-1}, {1}, {0}, {0;0.7854;1.5708}, {1,0.7071,0;0,-0.7071,-1})

Expected output:

Result
0	1
0.3927	0.9239
0.7854	0.7071
1.1781	0.3827
1.5708	0

Example 2: Scaled harmonic oscillator with sine solution

Inputs:

y_double_prime_coeffs	bc_ya	bc_yb	x_mesh	y_init
-4	0	1	0	0	0.7071	1
			0.3927	2	1.4142	0
			0.7854

Excel formula:

=SOLVE_BVP({-4}, {0}, {1}, {0;0.3927;0.7854}, {0,0.7071,1;2,1.4142,0})

Expected output:

Result
0	0
0.1963	0.3827
0.3927	0.7071
0.589	0.9239
0.7854	1

Example 3: Exponential growth with boundary conditions

Inputs:

y_double_prime_coeffs	bc_ya	bc_yb	x_mesh	y_init
0.5	1	2	0	1	1.2571	1.5843	2
			0.3333	0.5	0.5	0.5	0.5
			0.6667
			1

Excel formula:

=SOLVE_BVP({0.5}, {1}, {2}, {0;0.3333;0.6667;1}, {1,1.2571,1.5843,2;0.5,0.5,0.5,0.5})

Expected output:

Result
0	1
0.3333	1.2571
0.6667	1.5843
1	2

Example 4: Two coupled harmonic oscillators with boundary conditions

Inputs:

y_double_prime_coeffs		bc_ya	bc_yb	x_mesh	y_init
-1	0	1	0	0	1	0.7349	0.3888	0
0	-1	0	1	0.3333	-1	-1	-1	-1
				0.6667	0	0.3888	0.7349	1
				1	1	1	1	1

Excel formula:

=SOLVE_BVP({-1,0;0,-1}, {1;0}, {0;1}, {0;0.3333;0.6667;1}, {1,0.7349,0.3888,0;-1,-1,-1,-1;0,0.3888,0.7349,1;1,1,1,1})

Expected output:

Result
0	1	0
0.3333	0.7349	0.3888
0.6667	0.3888	0.7349
1	0	1

Python Code

import numpy as np
from scipy.integrate import solve_bvp as scipy_solve_bvp

def solve_bvp(y_double_prime_coeffs, bc_ya, bc_yb, x_mesh, y_init):
    """
    Solve a boundary value problem for a second-order system of ODEs.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        y_double_prime_coeffs (list[list]): Square coefficient matrix A in the equation y'' = A @ y (must be n×n)
        bc_ya (list[list]): Boundary condition values at the left endpoint x[0] (column vector with n rows)
        bc_yb (list[list]): Boundary condition values at the right endpoint x[-1] (column vector with n rows)
        x_mesh (list[list]): Initial mesh points for the solver (column vector, strictly increasing, minimum 2 points)
        y_init (list[list]): Initial guess for the solution (2n rows × m columns, where m is the number of mesh points)

    Returns:
        list[list]: 2D list [x, y1, y2, ...], or error message string.
    """
    def _ensure_matrix(data, name):
        """Convert input to a 2D matrix, validating structure and numeric content."""
        if isinstance(data, (int, float)):
            if not np.isfinite(data):
                return f"Invalid input: {name} contains non-finite values."
            return [[float(data)]]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        matrix = []
        row_length = None
        for i, row in enumerate(data):
            if not isinstance(row, list) or len(row) == 0:
                return f"Invalid input: each row in {name} must be a non-empty list."
            try:
                float_row = [float(value) for value in row]
                if not all(np.isfinite(value) for value in float_row):
                    return f"Invalid input: {name} contains non-finite values in row {i + 1}."
            except (TypeError, ValueError) as e:
                return f"Invalid input: {name} must contain numeric values (error in row {i + 1})."
            if row_length is None:
                row_length = len(float_row)
            elif len(float_row) != row_length:
                return f"Invalid input: all rows in {name} must have the same length (row {i + 1} has {len(float_row)}, expected {row_length})."
            matrix.append(float_row)
        return matrix

    def _ensure_column(data, name):
        """Convert input to a column vector (1D array), validating structure and numeric content."""
        if isinstance(data, (int, float)):
            if not np.isfinite(data):
                return f"Invalid input: {name} contains non-finite values."
            return [float(data)]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        column = []
        for i, row in enumerate(data):
            if not isinstance(row, list) or len(row) != 1:
                return f"Invalid input: each row in {name} must contain exactly one value (row {i + 1} has {len(row) if isinstance(row, list) else 'invalid type'})."
            try:
                value = float(row[0])
                if not np.isfinite(value):
                    return f"Invalid input: {name} contains non-finite values in row {i + 1}."
                column.append(value)
            except (TypeError, ValueError) as e:
                return f"Invalid input: {name} must contain numeric values (error in row {i + 1})."
        return column

    # Validate coefficient matrix
    matrix_result = _ensure_matrix(y_double_prime_coeffs, "y_double_prime_coeffs")
    if isinstance(matrix_result, str):
        return matrix_result
    coeff_matrix = np.array(matrix_result, dtype=float)
    if coeff_matrix.ndim != 2:
        return "Invalid input: y_double_prime_coeffs must be a 2D list (matrix)."
    if coeff_matrix.shape[0] != coeff_matrix.shape[1]:
        return f"Invalid input: y_double_prime_coeffs must be a square matrix (got {coeff_matrix.shape[0]}×{coeff_matrix.shape[1]})."

    # Validate boundary conditions
    bc_ya_result = _ensure_column(bc_ya, "bc_ya")
    if isinstance(bc_ya_result, str):
        return bc_ya_result
    bc_a = np.array(bc_ya_result, dtype=float)

    bc_yb_result = _ensure_column(bc_yb, "bc_yb")
    if isinstance(bc_yb_result, str):
        return bc_yb_result
    bc_b = np.array(bc_yb_result, dtype=float)

    n_equations = coeff_matrix.shape[0]
    if bc_a.size != n_equations:
        return f"Invalid input: bc_ya must have {n_equations} rows to match y_double_prime_coeffs dimension (got {bc_a.size})."
    if bc_b.size != n_equations:
        return f"Invalid input: bc_yb must have {n_equations} rows to match y_double_prime_coeffs dimension (got {bc_b.size})."

    # Validate mesh points
    x_mesh_result = _ensure_column(x_mesh, "x_mesh")
    if isinstance(x_mesh_result, str):
        return x_mesh_result
    if len(x_mesh_result) < 2:
        return "Invalid input: x_mesh must contain at least 2 points."
    x = np.array(x_mesh_result, dtype=float)
    x_diffs = np.diff(x)
    if not np.all(x_diffs > 0):
        # Find first non-increasing point
        bad_idx = np.where(x_diffs <= 0)[0][0]
        return f"Invalid input: x_mesh must be strictly increasing (x[{bad_idx + 1}] = {x[bad_idx + 1]:.6g} is not greater than x[{bad_idx}] = {x[bad_idx]:.6g})."

    # Validate initial guess
    matrix_result = _ensure_matrix(y_init, "y_init")
    if isinstance(matrix_result, str):
        return matrix_result
    y = np.array(matrix_result, dtype=float)

    expected_rows = 2 * n_equations
    if y.shape[0] != expected_rows:
        return f"Invalid input: y_init must have {expected_rows} rows (2 × {n_equations}) to match the expanded system (got {y.shape[0]})."
    if y.shape[1] != x.size:
        return f"Invalid input: y_init must have {x.size} columns to match x_mesh size (got {y.shape[1]})."

    # Define the right-hand side function for the first-order system
    # The second-order system y'' = A @ y is converted to:
    # u' = [y'; A @ y] where u = [y; y']
    def rhs(x_val, y_val):
        """Right-hand side of the first-order ODE system."""
        n_comp = y_val.shape[0] // 2
        result = np.zeros_like(y_val)
        result[:n_comp, :] = y_val[n_comp:, :]  # dy/dx = y'
        result[n_comp:, :] = coeff_matrix @ y_val[:n_comp, :]  # dy'/dx = A @ y
        return result

    def bc_residuals(ya, yb):
        """Boundary condition residuals at both endpoints."""
        residuals = np.concatenate([ya[:n_equations] - bc_a, yb[:n_equations] - bc_b])
        return residuals

    # Solve the boundary value problem
    # tol=0.001 provides reasonable accuracy for most engineering applications
    # max_nodes=1000 prevents excessive memory usage while allowing mesh refinement
    try:
        result = scipy_solve_bvp(rhs, bc_residuals, x, y, tol=0.001, max_nodes=1000)
    except Exception as exc:
        return f"Solver error: {str(exc)}"

    if not result.success:
        return f"Solver error: {result.message}"

    # Extract solution (only the first n components, not the derivatives)
    x_out = result.x
    y_out = result.y[:n_equations, :]

    # Validate that the solution is finite
    if not np.all(np.isfinite(x_out)):
        return "Solver error: solution x-coordinates contain non-finite values."
    if not np.all(np.isfinite(y_out)):
        return "Solver error: solution y-values contain non-finite values."

    # Format output as [x, y1, y2, ..., yn] for each mesh point
    combined = np.column_stack((x_out, y_out.T))
    return combined.tolist()

Overview

Excel Usage

Examples

Python Code

Online Calculator